#ifndef am_QsiCBW_H
#define am_QsiCBW_H

#include <string>
#include <math.h>
#include <tbb\concurrent_vector.h>
#include <tbb\blocked_range2d.h>
#include <boost\date_time.hpp>

#include <borealis\Raster.hpp>
#include "..\AuroraTypes.hpp"


namespace Aurora
{
	/*
	Class: am_QsiCBW
		Calculates the daily average incoming shortwave radiation following Walter 2005.
		CBW = Campbell, Bristow, Walter
		Latitude calculations are currently HARDCODED to Marmot Creek, Alberta. This WILL BE CHANGED IN THE FUTURE.
		The equation is based upon the solstices. Using the 2009 Solstice value of the 21st
	In Variables:
		Borealis::Raster<float> Girdded DEM
		Borealis::Raster<float> Gridded Max-Min Temperature (run am_TempInterp on Dt)
		string		  Date header name
	Out Variables:
		Borealis::Raster<float> Qsi
	References:
		Walter, M. Todd, Erin S. Brooks, Donal K. McCool, Larry G. King, Myron Molnau, and Jan Boll. 2005. Process-based snowmelt modeling: does it require more input data than temperature-index modeling? Journal of Hydrology 300: 65-75.
	*/
	class am_QsiCBW : public AuroraTypes::ModuleBase
	{
	private:

		//in
		int DEM;
		int GRID_DT;
		int DATE;

			
		//out
		int QSI;


	public:

		am_QsiCBW(std::string n)
		{
			DEM			= 0;
			GRID_DT		= 1;
			DATE		= 2;

			QSI			= 0;

			name = n;
		}

		class CalcQsi
		{
			 float const m_B;
			 float const m_So;
			 float const m_So30;
			 boost::shared_ptr<Borealis::Raster<float> > const m_Dt;
			 boost::shared_ptr<Borealis::Raster<float> > const m_Qsi;


		public:
			CalcQsi(float B, float So, float So30, boost::shared_ptr<Borealis::Raster<float> > Dt, boost::shared_ptr<Borealis::Raster<float> > Qsi)
				: m_B(B), m_So(So), m_So30(So30), m_Dt(Dt), m_Qsi(Qsi)
			{

			}

			void operator()(const tbb::blocked_range2d<unsigned int>& r) const
			{

				float B = m_B;
				float So =  m_So;
				float So30 = m_So30;
				boost::shared_ptr<Borealis::Raster<float> > Dt = m_Dt;
				boost::shared_ptr<Borealis::Raster<float> > Qsi = m_Qsi;
				


				for(unsigned int x = r.rows().begin(); x!=r.rows().end();++x)
				{
					for(unsigned int y = r.cols().begin(); y!=r.cols().end(); ++y)
					{
						float dt = Dt->GetRasterAt(x,y);

						//transmittance
						float expresult = exp( (-B/So30) * pow(dt,2.0F) );
						float Tt = 0.75F * (1-expresult);
						float pointQsi = So/1000.0F * Tt; //need to convert to MJ
						Qsi->SetRasterAt(x,y,pointQsi*1000000/(3600*24)); //conversion from MJ/m^2 to W/m^2
					}
				}
			}
		};

		void run(AuroraTypes::StationList* stations, AuroraTypes::VariableList* variables, AuroraTypes::InVariables& inVars, AuroraTypes::OutVariables& outVars)
		{
		
			//pi constant
			float PI = 3.14159F;

			//latitude in radians
			//50 N
			float phi = 0.889364355F;

			//day number
			boost::gregorian::date date = stations->at(0)->Now().GetGregorian();
			int J = date.day_of_year();
			
			//solar declination in radians

			float delta = 23.5F*PI/180*sin(((2*PI)/365)*(J-80));
			float delta30 = 23.5F*PI/180*sin(((2*PI)/365)*(J-110));//extra -30 is to get 30days previous for use below

			//117.5*10^3 kJ/(m^2 day)
			float Sprime = 117500;

			//max potential extraterrestrial solar radiation
			float t1 = acos(-tan(delta)*tan(phi)) * sin(phi) * sin(delta);
			float t2 = cos(phi) * cos(delta) * sin(acos(-tan(delta)*tan(phi)));
			float So =  (Sprime/PI)*(t1+t2 );

			float t130 = acos(-tan(delta30)*tan(phi)) * sin(phi) * sin(delta30);
			float t230 = cos(phi) * cos(delta30) * sin(acos(-tan(delta30)*tan(phi)));
			float So30 = (Sprime/PI)*(t130 + t230);
			So30 /= pow(10.0F,3.0F);
			
			boost::shared_ptr<Borealis::Raster<float> > Dt;
			boost::shared_ptr<Borealis::Raster<float> > Qsi;
			boost::shared_ptr<Borealis::Raster<float> > Dem;


			try
			{
				Dt = GetRaster(variables, inVars[GRID_DT])	;
				Dem = GetRaster(variables,inVars[DEM]);


				Qsi.reset(new Borealis::AsciiRaster(Dem->GetRow(),Dem->GetCol()));
				Qsi->SetXllCorner(Dem->GetXllCorner());
				Qsi->SetYllCorner(Dem->GetYllCorner());
				Qsi->SetCellSize(Dem->GetCellSize());
			}
			catch (std::runtime_error e)
			{
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + e.what());
			}
			

			boost::gregorian::days days(90);

			//these are the 2009 values
			boost::gregorian::date summerSolstice(2009,boost::date_time::Jun,21);

			boost::gregorian::date strtOfSummer = summerSolstice - days;
			boost::gregorian::date endOfSummer = summerSolstice + days;

			//assume it's summer
			bool isSummer = true;
			
			//if we lie outside of summer solstice +- 90 days, it's "winter"
			if(date < strtOfSummer || date > endOfSummer)
				isSummer = false;

			float B = 0.0F;
			if(isSummer)
				B = 0.282F * pow(phi,-0.431F);
			else
				B = 0.170F * pow(phi,-0.979F);

			tbb::blocked_range2d<unsigned int> range(0,Dem->GetRow(),0,Dem->GetCol());
			tbb::parallel_for(range,CalcQsi(B,So,So30,Dt,Qsi));

			AuroraTypes::VariableList::accessor a;
			variables->insert(a,outVars[QSI]);		
			a->second = Qsi;

		}



	};

}

#endif